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Abstract. - An investigation is made of multihump and periodic solutions of the semiclassi- 
cal coupled equations describing laser radiation copropagating with a Bose-Einstein condensate. 
Solutions reminiscent of optical vector solitons have been found and have been used to gain un- 
derstanding of the dynamics observed in the numerical simulations, in particular to shed light on 
the phenomenon of jet emission from a condensate interacting with a laser. 
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Introduction. — The experimental realisation of 
Bose-Einstein condensates (BEC) and of coherent atomic 
beams has inspired a renewed effort in the theoretical mod- 
elling of such physical systems. Of particular interest from 
the point of view of applications, is the possibility of ma- 
nipulating atomic beam structures by using their interac- 
tions with light. Besides, due to the mathematical mod- 
elling of these interactions, these phenomena can also offer 
a possible test of the analogies between optics and quan- 
tum matter waves, which is what we would like to explore 
further in this work. 

It is well known that the dynamics of a BEC satisfies 
an equation which has the same mathematical structure 
of the nonlinear Schrodinger equation, the fundamental 
equation of nonlinear optics, [T], the Kerr-like nonlinear- 
ity being due to the atom-atom interactions which can 
be either focusing or defocusing depending on the sign of 
the scattering length. Considering the interactions of a 
BEC with laser radiation brings new terms into the atom 
equation whose form depends on the modelling of the in- 
teractions, the simplest form being the dipole-dipole one. 
Saffman started to study the consequences of these inter- 
actions, [2] while Krutistky et al. derived a fundamental 
model from first principles within the framework of quan- 
tum field theory and later generalized it to consider atom 
transitions (in the case of high atom densities) , [3] ■ More 
recently, the same structure for the atom equation was 
derived within a semi-classical approach, [4]. The idea 



there was to work out an expression for the force exerted 
by light on atoms from basic classical physics and use it 
as a potential term in the Schrodinger equation of the 
atoms. It is important to notice that, within the limits of 
a semi-classical approach, the atom equation agrees with 
the more general equation of [3] . The point is that a natu- 
ral consequence of BEC atoms interacting with photons is 
the emergence of a coupling between the dynamics of the 
atoms and that of the laser radiation. In fact, what these 
models describe is photons from the incident radiation ex- 
citing atoms which in turn re-emit photons, which are 
then absorbed by other atoms, thus giving rise to a long- 
range interatomic interaction, [HIE]. This obviously has 
an effect not only on the atoms but also on the radiation. 
This gives rise to a coupling term in the atom equation 
and the need for an additional equation describing also 
the dynamics of the radiation. It is this coupling that can 
bring about novel effects. It has been shown via numeri- 
cal simulation of the coupled dynamics that self-localised 
structures and mutual atom-light guiding can be achieved 
exploiting the dipole-dipole interactions, [6]. In particular 
it was observed that solitary-like localised structures can 
be generated and emitted from a central bunch of atoms, 
[7] . This is suggestive of analogous optical effects such as 
soliton ejection, see [5] and references therein, and we will 
try here to obtain an understanding of how the structures 
ejected in laser-BEC interactions are created and emitted. 
The peculiarity of the emission discussed here is that the 
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generation of the solitary-like structures is entirely due 
to the atom-laser coupling and their ejection from the re- 
gion where they have been created can be ascribed to their 
tail-tail interaction. This sort of interaction, being the 
structures initially well separated, is determined by the 
overlap of the structures themselves, i.e. the cross terms 
arising when calculating terms such as \iPtot\ 2 if iPtot 
can be written approximately as the sum of two or more 
localized functions. Since these functions are assumed to 
be spatially separated, the only contributions to the in- 
teraction term |?/>tot| 2 will come from their overlapping 
tails. 

In optics, for example in the early study of Gubbels et al., 
[9], up to the more recent work of Assanto et al., [ID] , the 
emission is engineered in such a way that the trap or guide 
does not undergo any modification. Similarly in studies of 
soliton emission from BECs, see for example [TT], the ba- 
sic physical mechanism for emission concerns the atom- 
atom interactions and the trap is not modified. In our 
case, the atoms modify the medium through which they 
propagate, i.e. the laser radiation, which in turn affects 
the atom propagation via a nonlinear coupling. If mutual 
trapping can be achieved, localised symbiotic atom-light 
structures will start to propagate together and the emis- 
sion will be due to the interaction of these newly born 
structures rather than to tunneling or engineering of the 
trap. This process is also reminiscent of the solitonic glu- 
ons studied by Ostrovskaya et al., [12] . 
Inspired by numerical simulations which clearly show how, 
during the coupled propagation of laser and atoms, sym- 
metric and mutually localised structures are formed and 
then ejected, the first step we would like to take is an inves- 
tigation into the possible existence of multi-hump station- 
ary solutions of the coupled equations. This kind of solu- 
tion is not found in the usual nonlinear Schrodinger equa- 
tion but is known to be a possibility in coupled systems of 
nonlinear equations such as those describing birefringent 
fibers. If multi-hump solutions exist for the laser-BEC 
system as well, it should be possible to use known methods 
of nonlinear optics and to study the interaction between 
the different soliton-like formations to infer whether there 
will be mutual repulsion or attraction or whether they can 
propagate together. This would justify the numerical re- 
sults and provide an insight into the possibility of using 
these effects to manipulate localised atom-light structures. 
In what follows we will briefly review the basic coupled 
model and stationary coupled equations as well as the 
limitations imposed by a semi-classical description. We 
will then search for multi-hump solutions and study in 
particular the tail interaction of double peaked solutions 
to show how, within this model, it is possible to expel 
solitary-like structures in a jet-like fashon, or even having 
them colliding and merging together into a single central 
structure. 



Model equations. — The semi-classical coupled sta- 
tionary equations as introduced in [HE] are given by 
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where for the full atom wave function and the full laser 
field it has been assumed 



¥(r,i) = $(r)exp(-twot), 
E(r,i) =Re[S(r)exp(-iu} L t)]. 



(3) 
(4) 
(5) 



Here, a(u>) = —d 2 /TiA is the atomic polarizability at the 
laser frequency ui = k^c, with A = ujl — uj a being the 
detuning from the nearest atomic resonance frequency oj a , 
and d is the dipole matrix element of the resonant tran- 
sition. Besides, for the gas density we have n — \^\ 2 , 
Uq — 4nh 2 a s /m, m is the atom mass and a s is the s— 
wave scattering length (which will be assumed positive as 
for repulsive atom-atom interactions) and Hq is the linear 
single-particle Schrodinger Hamiltonian. Furthermore, in 
order to reduce Maxwell's equations to three scalar equa- 
tions, it has been assumed that L n ^> Al and Ve • E ~ 
(L n is the characteristic length scale of transverse density 
modulations and Al is the radiation wavelength). Notice 
that the atom equation is the approximate version of what 
was presented in [3], valid only under the assumptions of 
a semi-classical approach which limits the model to a well 
defined range of parameters: The concept of force being 
purely classical, quantum fluctuations, stochastic heating 
and any incoherent process are to be neglected which is 
a consistent assumption if large detunings are considered 
|A| 3> uj a , r (r is the natural line width of the atoms). Al- 
though these limitations are quite strict, the model should 
elucidate the basic physics of the interaction and hopefully 
the resulting structures will be resilient enough to be in- 
teresting even under non ideal conditions. 
Finally, for a mean field model to be valid for the atom 
wave function, we must consider not only the zero tem- 
perature limit but also a low density limit with na^ <C 1, 
see pp. Furthermore, a low density regime is required in 
order to avoid the singularity of the model and concomi- 
tant spurious collapse-like phenomena. In view of this 
limitation and in order to simplify the analysis, in what 
follows we will neglect the denominator in both equations. 
This is a delicate step and cannot be taken when analysing 
the coupled dynamics of the system, since it is not guar- 
anteed that during the evolution the peak values of the 
atom density will satisfy the low density assumption (as 
seen in numerical simulations of the coupled evolution, [4] , 
the nonlinear focusing action of the laser could be strong 
enough as to focus the atoms to very high peak densities 
in a single spot thus breaking the low density assump- 
tion). However, it is possible to accept this simplification 
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when studying the stationary solutions with the proviso 
that only low density solutions will be accepted. 

The physical setup of the problem. In general, 
if laser and atoms are prepared in an initial state which is 
not a fully stationary state of the system, we expect the 
system to evolve showing changes in both the atom wave 
function and the laser amplitude profile with the propa- 
gation variable (hereafter chosen to be z).Thus it was ob- 
served in [6] and [7] via numerical simulations of the cou- 
pled equations that, depending on the parameters of the 
initial state, the propagation effects could have different 
outcomes. In particular, starting from a super-Gaussian 
laser amplitude much wider than the initial Gaussian atom 
density distribution, a regime of parameters (peak laser in- 
tensity and peak atom density) could be found for which 
the initial single-peak atom wave function and the flat 
laser amplitude profile slowly changed into two symmet- 
rical peaks. Once these structures were formed, further 
propagation led to different possible scenarios: (1) out- 
ward motion of the two peaks moving farther apart from 
each other (this is what we call "jet emission"); (2) inward 
motion of the two peaks resulting in coalescence into a sin- 
gle central peak; or (3) inward motion of the two peaks 
resulting into a bound state with the two structures oscil- 
lating about the central position. Thus, from simulations 
of the coupled propagation equations it seems possible to 
see the creation of such symmetrical double-hump struc- 
tures which however do not survive for long and undergo 
their own motion. In order to understand their nature and 
destiny (can they be thought of as solitary waves? Can 
we describe the outcomes of their coupled propagation in 
terms of soliton-soliton interactions?), we need to know 
what sort of stationary structures the model admits and 
work out the properties of their interactions. 

Stationary solutions. — The existence of stationary 
solutions and their stability is fundamental from the point 
of view of realisable structures. We are interested in sta- 
tionary solutions corresponding to mutual guiding in the 
form of long distance localised beam propagation. To anal- 
yse these solutions we focus on modes that are localised 
in the transverse direction, the analogue of Kerr spatial 
solitons, by assuming 

£{r) = a(r±) exp(ih a z)e (6) 
$(r) = <f>(r±)exp(ihfz) (7) 

where rj_ denotes the dimension transverse to the propaga- 
tion direction z, e is the polarization vector of the electric 
field, a, 4> are real amplitudes and h a , are the laser and 
atom propagation constants respectively. These solutions 
must satisfy the fully stationary equations derived from 
([1]) and @, which when written in normalised variables 
become: 

V 2 ± 4> - f34> 3 + sa 2 </> - 2^k^ = 0, (8) 
Via + 3s(f) 2 a - 2n a a = 0. (9) 



The normalisation used is: f = rfct, for the atom wave 
function <fr = 4>/<fi* with (47r|a|/3)0j = 1, for the laser 
a — a/a*, with m\a\a 2 / {2Ti 2 k\) = 1, s = sign(a), 
and (3 — 6a s /(kj j \a\) is directly proportional to the 
strength of the collisional nonlinearity. The tilde will be 
dropped hereafter unless otherwise stated. Furthermore, 
/j = k a /k L , k a = y / 2muj a /h, kj, = (ti^/kl - ij /2 and 

K a = (b? a jk 2 L — l) /2. For simplicity, we will consider only 
one transverse dimension, r_L = x and assume fj, = 1. To 
allow for mutual trapping, we will also assume s = +1. 
A further change of variables allows to rewrite these equa- 
tions in terms of the relative wave number only 

4>" + (a 2 - 1)0 - /30 3 /3 = 0, (10) 
a" + {4> 2 - R)a = 0. (11) 

where R — K a j k§ is the relative wave number, x = 
x^2fi 2 K rt „ a 2 = a 2 j (2^ 2 ^), cj) 2 = 30 2 / (2fj 2 n l p). Again 
the bar will be dropped hereafter. 

To find localised solutions, we can solve this set of equa- 
tions numerically as an eigenvalue problem. A first inte- 
gral of motion of (fTU)) and (fTTj) can be obtained by mul- 
tiplying (fTU)) by (j)' and (ITT1) by a' and integrating over x. 
Combining the two resulting integrals we find 

4>' 2 + a' 2 ~ <t> 2 - na 2 + cf> 2 a 2 - -f3<p 4 = constant, (12) 

6 

The constant is zero for a localized solution vanishing at 
infinity. From this, we can deduce a relation between the 
values of </> and a at the symmetry point x — where we 
require 0'(O) = a'(0) = 0: 




where ao = a(x = 0),4>o = 4>(x = 0). Then ao can be 
used as a shooting parameter to search for localised so- 
lutions with the shooting method [13]. For a fixed value 
of the eigenvalue n, a value is assumed for ao with <ji>o 
consequently calculated from (fT"2")) and a solution is calcu- 
lated with an ordinary differential equation solver. The 
value of ao is than varied until a solution with the de- 
sired characteristics is found (single-hump, double-hump 
...). As shown in our previous work [1], this equation 
admits single hump soliton-like solutions, in particular, 
it was shown that for one-scaled distributions (n = 1) 
the set of eqs. (|10[) . (fTTj) reduces to one equation if 
a{x) = \J /3/3 + 10(x), the soliton solution correspond- 
ing to ao = V2 as is evident when solving the resulting 
nonlinear Schrodinger equation. We have expanded the 
family of solutions of (fTU)) . (JTTJ) and it is now clear that 
this set of equations admits multi-hump and periodic so- 
lutions as well for which the Hamiltonian (fT2|) is exactly 
equal to zero. This can be easily seen in the simple case 
of j3 = 0. At k ^ 1 with a = \/2, the coupled equations 
admit periodic solutions, see for instance fig[T] 
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Fig. 1: Periodic solutions for ft — at different values of k. 
Blue line 4>(x), red line (always positive) a(x). All quantities 
normalized as in the text. For re > 1 the role of laser and atoms 
is simply inverted. Here k$ = 10 -3 . Colour online. 



Multi-hump solutions can be found as well, fixing re 
and varying aQ. Some examples for /3 = are shown in 
figl2ja) and (b) while part (c) and (d) show examples of 
different multi-hump solutions for /3 ~ 38 corresponding 
for instance to a detuning of 100 times the decay rate for 
87 Rb atoms (the same value used in the previous numerical 
simulations). By varying the shooting parameter do it is 
possible to find solutions with different number of peaks, 
thus it seems indeed possible to find multi-hump solutions 
for any value of /?. 




Fig. 2: Multi-hump solutions at re = 1.01, /3 = [(a) and (b)] 
and 38.4281 [(c) and (d)]. Blue line (always positive in (a) 
and (b)) <f>(x), red line (always positive in (c) and (d) a(x). 
All quantities normalized as in the text. Shooting on the peak 
laser amplitude ao it is possible to identify different solutions 
with increasing number of peaks. Here = 10 -3 . Colour 
online. 



Jet emission. — We would like to underline the fact 
that these kind of multi-peak structure was never clearly 
seen in numerical simulations of the coupled propagation 
equation. As mentioned in the sectiom "The physical 
setup of the problem" , we mainly observed the generation 
of two symmetric peaks (or four at higher initial atom 
densities) thus we cannot say anything about the gener- 
ation of many-peak solutions within the range of param- 
eters we have explored (the main limit being low atom 
densities). It is very likely that they can be generated 
at higher atom densities but appear only as a transient 
to transform quickly into more stable structures with two 
peaks only. This has obvious consequences for what could 
be oserved experimentally therefore it is important to un- 
ravel the dynamics of the structures. In particular we shall 
study the simplest case of two symmetric structures as the 
model case for this kind of dynamics. A knowledge of how 
they evolve gives the only acceptable indication of what 
could be observed. Therefore, once the existence of multi- 
hump solutions has been established, we can proceed to 
study the existence of stationary solutions that originate 
the ejected jets. The simplest case found in numerical sim- 
ulations, [7], pictures an initially centrally localised bunch 
of atoms coupled to a localised profile of laser intensity (a 
Gaussian and a super-Gaussian respectively in the simula- 
tions), splitting into two "jets" which are then ejected and 
start to propagate in opposite directions together with two 
laser jets. The aim now is to find stationary solutions con- 
sisting of only two symmetric (or antisymmetric) localised 
humps. If they exist, we expect a tail-tail interaction be- 
tween them because of the nature of the equations. In 



the case of well separated humps, i.e. weak overlap and 
consequently weak tail interaction, a perturbative analysis 
completely akin to that introduced by Yang in [14] can be 
used to make predictions about the fate of the two humps. 
The question is whether they will propagate together or 
whether they will repel or attract each other because of 
their interaction. 

An example of a symmetric two-peak solution is shown 
in figE] for fixed /3 and re. These structures are found 
by imposing cf>(0) = a'(0) = 0, using the Hamiltonian 
(fl~2"]) to obtain an expression for a(0) 2 = 0'(O) 2 /re and 
shooting on the value of the first derivative 4>'(Q). Due 
to the symmetry of the equation, it is enough to cal- 
culate one half of the solution, say for x > and any 
combination (laser symmetric-atoms antisymmetric, laser 
antisymmetric-atoms symmetric and so on) is still a solu- 
tion. 

With a known double-peak solution showing two very 
well separated peaks, it is possible to study the nature of 
the interaction between the tails of the two peaks. We 
have followed the procedure by Yang in |14) . who intro- 
duced a perturbative analysis of the interaction of vector 
solitons which we shall only briefly recall here. Modify- 
ing the perturbative approach of Karpman and Solov'ev, 
Yang found dynamical equations for the parameters of 
two-vector solitons weakly interacting with each other be- 
cause widely separated. In particular, for cases in which 
the two components have phase difference or n, Yang 
determined an equation for the dynamics of the vector- 
soliton separation Aa; in the form of a Newtonian equa- 
tion of motion md 2 Ax/dz 2 = — W. The fixed points of 
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Fig. 3: An example of a double-hump solutions at k = 1.1 with 
P — 38.4281. Blue line <j>(x) (antisymmetric curve at lower am- 
plitude), red line a(x) (symmetric curve at higher amplitude). 
All quantities normalized as in the text. Here = 1CP 3 . 
Colour online. 



this equation give the separation distance of the two com- 
ponents of a stationary solution. We can therefore apply 
Yang's formulation to the main case under analysis, that 
of a two-peak stationary solution as the one numerically 
evaluated and shown in figEl This solution corresponds 
to what seen in numerical simulations of the atom-laser 
coupled propagation (see the section "The physical setup 
of the problem" ) before the two peaks start to move either 
inward or outward. With a detuning of about 100 times 
the decay rate for 87 Rb atoms (i.e., (3 ~ 38), we find a 
peak atom density of 2.49 10 19 m -3 and a peak laser in- 
tensity of 0.051 mW/cm 2 . Utilizing Yang's formula, we 
arrive for the case shown in fig|3] at Aio = 1.3344 10 3 
which compares quite well with the separation found nu- 
merically 1.3359 10 3 . Therefore we can use the potential 
V to infer the dynamics of the separation Ax. Again for 
the case shown in fig(3]and using Yang's formula for V, the 
potential is shown in figlU It is clear that if the two com- 
ponents are generated at a distance Ax < Axo, the two 
solitons will repel each other and be ejected away whereas 
two components generated at Ax > Aio will attract each 
other with the possible creation of a bound state where 
the two components oscillate about their equilibrium po- 
sition. This model, based on the existence of multi-peaks 
solutions of our stationary equations, thus explains quite 
well the qualitative features of the numerical results of [7] , 
both the jet emission and the formation of a bound state. 
It thus seems possible to say that the dipole-dipole inter- 
action can lead to the generation of two-peak structures 
for both atoms and laser but the destiny of the two peaks 
is that of moving apart from each other in a solitary-like 
fashion or to move towards each other. This is indeed 
what seen via numerical simulations of the coupled prop- 
agation and we interpret the two-peak solutions as the 
unstable seeds of the evolution predicted by Yang's model 
and observed numerically. 

Conclusions. — Numerical simulations of copropa- 
gating laser-BEC systems seem to indicate that not only 



Fig. 4: Yang's potential V versus the soliton separation Ax for 
k = 1.1 and = 38.4281. 



mutually localised structures can be formed, but also that 
solitary-like wave packets can be emitted out of the cen- 
tral interaction region. In order to understand the physics 
of the jet ejection, we have used a simplified system of 
coupled equations and studied the existence of possible 
multi-peaked stationary solutions, usually not possible for 
the nonlinear Schrodinger equation. Applying the results 
of a perturbative method elaborated in nonlinear optics to 
study their dynamics, we have found that jet emission and 
bound state formation can be explained by the tail inter- 
action of the soliton components of a two-peak structure 
within the model we are using. 
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